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1 INTRODUCTION 

The objective of this numericalVexperimental study is 
to improve the understanding of the effects of gravity 
on the two-way interaction between dispersed particles 
(bubbles or liquid droplets) and the carrier turbulent flow. 
Due to the imposed page limit, only a condensed de- 
scription of the results will be presented. Reference is 
made to other publications for the details. 

The experiments reported here are concerned with 
the dispersion of liquid droplets by homogeneous tur- 
bulence under various gravitational conditions and the 
effect of these droplets on the evolution of the turbu- 
lence of the carrier fluid (air). 

The numerical results presented here are obtained 
from direct numerical simulations (DNS) of bubble - 
laden turbulent flow to examine the effects of small bub- 
bles on the behavior of decaying turbulence. All pub- 
lished DNS studies of turbulent flows laden with parti- 
cles, except [7], adopt the Lagrangian-Eulerian approach. 
In this approach the carrier flow velocity field is ob- 
tained by solving the Navicr-Stokes and continuity equa- 
tions at fixed mesh points whereas the trajectories of 
the dispersed particles are computed by solving the La- 
grangian equation of particle motion f 1 ],[2] . Using this 
method to simulate turbulent flows with two-way cou- 
pling between the particles and the carrier flow is lim- 
ited at present and in the foreseeable future by the mem- 
ory and speed of available supercomputers, including 
the latest parallel supercomputers. This limitation forces 
all DNS of particle - laden turbulent flows (with two- 
way coupling) to compute the trajectories of only a frac- 
tion of the actual number of particles. The accuracy of 
DNS results is directly proportional to the magnitude 
of this fraction, being highest when the fraction equals 
unity [21. 

The alternative approach for predicting particle-laden 
flows is known as the ‘two-fluid’ (TF), or ‘Eulerian- 
Eulerian’, approach [3], and has been used only with 
the Reynolds - averaged equations of motion, not with 
DNS. In the TF approach, the governing equations are 
obtained by volume - averaging the equations of motion 
of both phases (the carrier flow and particles) based on 
the assumption that the dispersed particles behave as a 


“continuum” under certain conditions. 

The objective of the numerical part of this paper is 
to describe briefly how DNS can be performed using 
the two-fluid approach for bubble - laden homogeneous 
isotropic turbulence without forcing. More details are 
given in [7]. Turbulent homogeneous shear flows laden 
with droplets/bubbles will be studied in the next phase 
of the project. 

2 NUMERICAL STUDY 

2.1 Governing Equations 

We consider spherical bubbles with diameter d much 
smaller than the characteristic length scale of the flow, 
Lj, and average the equations of motion of the fluid 
and bubble over a length scale scale A which is much 
smaller than L / but much larger than the bubble diame- 
ter, d < A < L j . Thus the bubble phase can be treated 
as a continuum characterized by the velocity I )(r,t) and 
concentration (or, volume fraction) C(r,l)= -r/ :i n(r,t)/6, 
where n(r,t) is the bubble number density. 

We assume that the density of the gas and, conse- 
quently, the mass of the bubble are negligible compared 
to those of the surrounding fluid, p; » pi = 0. We 
also assume that the bubble concentration, C, is small 
enough (i.e. C < 10~ 3 ) and thus neglect its contribu- 
tion to the fluid inertia and continuity, i.e. we retain C 
only in the buoyancy term in the momentum equation 
of the carrier flow. This is analogous to the Boussinesq 
approximation in a stratified fluid with effective density 

(1 - c ) P{ . 

The resulting equations of the conservation of the 
fluid and bubble phase momentum and mass are [7]: 

DU 1 

d t p+v± U i +(C-<C» gS i: , ( 1 ) 

Dt p } 

dj Uj = 0 , ( 2 ) 

^ = 3 ^ + -((/,— V,+WS, ; ), ( 3 ) 

at L)t Tb 

+ fyC Vj = 0 , ( 4 ) 
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where U, is the fluid velocity, and I is the velocity of 
the bubble phase. The Lagrangian derivatives D/Dt, — 
d/dt + Ujdj and d/dt = d/dt + Vjdj are taken along 
the trajectories of the fluid point and bubble, respec- 
tively, and g is the projection of the gravity acceleration 
on the z— axis, gi = —gS( : . The modified hydrostatic 
part P of the pressure field in eq.(l) is defined as: 

P - P + p f g [ (1 -<C»dz, (5) 

Jo 

where < C > is the ensemble-averaged bubble con- 
centration. We evaluate < C > as an average over a 
horizontal (z) plane. The bubble response time n and 
terminal velocity IT are defined as 



( 6 ) 


H’ = 2 r b g. (7) 


The momentum conservation and continuity equa- 
tions (l)-(4) for both phases are solved in a cubical do- 
main with periodic boundary conditions. The equations 
are discretized in an Eulerian framework using a second- 
order finite-difference technique on a staggered grid con- 
taining 96 3 points equispaced within unit length in each 
of three coordinate directions (x,y,z). The Adams - Bash- 
forth scheme is used to integrate the equations in time. 
Pressure is obtained by solving the Poisson equation us- 
ing Fast Fourier Transform. More details about the nu- 
merical method and its accuracy are discussed by EI- 
ghobashi and Truesdell [6], 


wavenumber, and i) is the Kolmogorov length scale rj = 
( 1 / 3 /e) 1/M . The bubble response time is restricted by the 
two conditions : 

dt < V > (9) 

Re b = — < 1, (10) 

v 

which are necessary for the derivation of (3) [7]. Substi- 
tuting the terminal velocity (7) and the bubble diameter 
(d b = (36 ^n) 1 / 2 ) in (9) and (10), and using the equal- 
ity if = utic, we rewrite the conditions (9) and (10) 
respectively as 


1 

< 36 “ 

0.028, 

(11) 


1/3 

(12) 

MW) 

= T t . 


Both conditions (11) and (12) are satisfied throughout 
our simulations. 

Figure 1 shows the time evolution of the concentra- 
tion spectrum. It is evident that no numerical instability 
occurs. The spectrum Ec{k) at high wave numbers ap- 
proaches an asymptotic form at t=10. The high wave- 
number range in the spectrum (i.e. k > 40) would de- 
tect any numerical instability if it existed. The reason 
for the absence of the instability is that the fluctuations 
of the bubble concentration, caused by the preferential 
accumulation, are proportional to the ratio which 
decreases with time (approximately as ~ 1/f) in decay- 
ing turbulence. 


2.2 Results 

In this section we present the DNS results for bubble 
dispersion in isotropic decaying turbulence with both 
one-way and two-way coupling. 

2.2.1 Dispersion of bubbles in isotropic decaying 
turbulence (with one-way coupling) 

DNS of bubble dispersion in isotropic decaying turbu- 
lence is performed with the initial condition: Re ao = 
25. The initial bubble velocity and concentration are 
prescribed as: 

I , — A, - TT , Co = Qo , (8) 

where the bubble terminal velocity IT is given by (7). 

The ability of the simulation to resolve the motion 
at the smallest turbulence scales is assured by the crite- 
rion ij k m ai > 1, where k ma x is the highest resolved 


Re X0 =25,x o =O.04T k0 (1-way,C 0 =a 0 ) 



Figure 1 : Time evolution of spectrum of bubble concen- 
tration fluctuation 
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It should be noted that both the accumulation of bub- 
bles and the absence of the diffusivity in the transport 
equation for the bubble concentration (4) may lead to 
instabilities in the numerical solution due to the devel- 
opment of steep gradients in the concentration field. The 
occurrence of this numerical instability depends on the 
initial distribution of the bubble concentration, the flow 
Reynolds number and the bubble response time. In our 
DNS we chose the initial microscale Reynolds number 
Re\o = 25, so that at time of the injection of bubbles 
(t= 1) the small-scale motions are resolved, i.e. k ma , T q > 
1 where k maT = N g n is the maximum wave number for 
the given grid resolution N g — 96. The numerical in- 
stability may occur for higher-inertia bubbles, i.e. for n 
of the order of the Kolmogorov time scale n- . However, 
prescribing rt, ~ 77. would violate the condition df, < t], 
which is necessary for deriving equation (3) of bubble 
motion. 

2.3 Two-way coupling effects on decaying 
turbulence 

Here we examine the effects of the dispersed bubbles 
on the temporal development of decaying isotropic tur- 
bulence. We consider three cases with different initial 
bubble concentration profiles in z-direction, but with the 
same bubble response time as in the one-way coupling 
case. 

The first case is for a uniform initial bubble concen- 
tration, C 0 = a 0 , where a 0 is a reference concentration 
set equal to 0.005 to allow neglecting bubblc-bubble in- 
teractions. The second case is for stable linear stratifica- 
tion, with a constant concentration gradient in the verti- 
cal (z) coordinate, C 0 = o 0 (l + e), while the third case 
is for unstable linear stratification, Co = «o(2 - z). In 
the cases of stable and unstable stratification, periodic 
boundary conditions in the z-direction are imposed on 
the instantaneous concentration fluctuation C' = C- < 
C > . We first consider the modification of the turbu- 
lence energy spectrum E(k), which is governed by 

0 t E{k) = T{k)-e(k) + 9 b (k), (13) 

where e (k) is the dissipation rate of E(k), and T(k) is 
the rate of energy transfer between the wave numbers 
(i.e. between the different scales). The source of the 
modification of the energy spectrum and spectral trans- 
fer process is $(,(&) which can be regarded as a spectral 
buoyancy flux, analogous to that in a stratified fluid with 
density C', and is defined as 

9 b {k) = gRe{C , (k)U:(k)}. (14) 


Re lo =25,T b =0.O4T ko ,t=3.0 (2-way) 



Figure 2: Modification of the turbulence energy spec- 
trum in a bubble-laden decaying turbulence 


Figure 2 shows the difference between the energy 
spectra in two-way and one-way coupling cases com- 
puted at lime t=3. As expected, the turbulence energy 
increases in the case of unstable stratification and is re- 
duced in the case of stable stratification. In the non- 
stratified two-way coupling case, the spectrum remains 
practically unchanged compared to the one-way cou- 
pling case. 

3 EXPERIMENTAL STUDY 

The objectives of the experiments reported here are mea- 
suring the dispersion of the droplets by homogeneous 
turbulence under various gravitational conditions and 
determining the effects of the presence of the droplets 
on the evolution of the turbulence of the carrier fluid 
(air). 

Experimental facility and initial conditions 

The experimental facility (Fig. 3) consists of a blow- 
down wind tunnel made up of two sections. The first 
section serves as the two-phase flow generation system, 
while the second is a 2 meter long test section of 20 x 
20 an square cross section. The air flow, which is sup- 
plied by a blower, passes through a slight contraction 
duct followed by a honeycomb and a grid which damp 
out the fluctuations generated by the blower. The air 
stream then flows through the spray generation system. 
This system consists of two arrays of 33 airblast atomiz- 
ers each arranged so that the air flow is uniformly laden 
with droplets. Each atomizer is made of a pressurized 
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Figure 3: Experimental facility. 


air jet impinging normally onto a water jet, thus cre- 
ating a spray whose main characteristics are given by 
the pressure of the air and the water flowrate supply- 
ing the liquid jet [ 1 0],[ 11]. In the experiments reported 
here, only one of the two arrays was supplied with wa- 
ter, thus only half of the cross section of the channel 
flow was loaded with droplets. A splitter plate is placed 
horizontally at the mid-plane of the channel, separating 
the droplet-laden flow from the non laden air flow (see 
Fig. 3). The length of the splitter plate is long enough 
to ensure that the size distribution of the droplets is uni- 
form across the cross section of the channel before en- 
tering the test section. Note that x, y. z correspond re- 
spectively to the longitudinal, normal and spanwise co- 
ordinates. Due to the spanwise uniformity of the flow, 
all measurements reported here were made at the z = 0 
plane. 

Preliminary measurements have been carried out to 
determine the characteristics of both the turbulent air 
flow and the spray. Velocity measurements were made 
using hot-wire anemometry. Measurements of mean 
longitudinal velocity profiles show that the small wake 
resulting from the presence of the splitter plate follows 
the common decay law ^ oc xr 1 / 2 , where MJ is 
the velocity deficit in the wake. Since for x > 35 cm, 
< 10%, most of the measurements reported here 
have been done beyond this downstream location where 
it is then considered that the wake has no longer influ- 
ence on the flow. As we are interested in the interac- 
tion between the particules and a homogeneous turbu- 
lent flow, we restricted our investigation to the region 
-40 mm < y < 40 mm in order to avoid any effect of 
the boundary layers forming at the top and bottom sur- 
faces of the test section. The mean velocity in the free 
stream selected for these experiments, U, is 14 rn/s, 
while the longitudinal turbulent intensity , u'/U is mea- 
sured to be 2%. This corresponds to mean and turbu- 
lent Reynolds numbers, Re = and Re, = i! ~, of 
1.9 10 s and 1900 respectively. H is the height of the 
channel, and L is the integral scale of the flow, calcu- 
lated from the turbulent spectrum measured in the free 



Figure 4: Probality density function for o. y = 

— 30; □, 0; o, 30 mm, — respective log-normal fits, 
x=19 cm, stable stratified case. 



y (mm) 


Figure 5: Flux of particles profiles, x=107 cm. o, D < 

5; □, 10; o, 15; x. 25; +, 35; A, 50; V, 65; •, 80 fim. 


stream. It is approximately 10 cm. The droplet void 
fraction in the free stream is a = 1.8 10 -5 . The spray 
size distribution and the velocity of the droplets were 
measured with a Phase Doppler Particle sizer (Aero- 
metrics PDPA). This instrument has already been tested 
successfully on a similar experiment [9]. The initial 
size distribution is shown in Fig. 4 by the open circles. 
The corresponding mean size is D io = 15 fim and the 
Sauter Mean Diameter is D 32 = 33 /im. As can be 
seen, the initial probability density function for the droplet 
sizes is well described by a log-normal distribution. As 
we are interested in determining the turbulent diffusivity 
of the droplets of various sizes, velocity measurements 
and concentration of droplets have been done for differ- 
ent sizes of particles. For this purpose, PDF have been 
discretized into 8 bins corresponding respectively to di- 
ameters, D, less than 5, 10, 15, 25, 35, 50, 65 and 80 
fim. 
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Figure 6: Mean and fluctuations velocity pro- 

files, x-67 cm. Open symbols, longitudinal ve- 
locity; closed symbols, normal velocity, o, D < 
5; □, 15; o, 35; A, 80 /rm 

3.1 Results 

Two sets of experiments were conducted to systemati- 
cally study the role of gravity. In all the experiments, 
only half of the wind tunnel entrance cross section was 
laden with uniform-concentration spray of droplets of 
a given polydispersed size distribution shown by the 
open circles in Fig. 4. In the first set, the uniform- 
concentration droplet spray initially occupied the lower 
half portion of the entrance cross section (stable strat- 
ified case), while in the second set the spray initially 
occupied the upper half portion of the wind tunnel cross 
section (unstable stratified case). 

Stably-stratified case 

The mean and fluctuations of the longitudinal and nor- 
mal velocities are shown in Fig. 6; for clarity, only four 
bins are shown. One observes that, at x = 67 cm there 
is still a slight influence of the splitter wake. Indeed, 
around the position y = 0 mm, there is a slight effect 
in the mean velocity and an increase in both fluctuating 
components of the velocity. Nevertheless, one can con- 
clude that the flow is almost homogeneous turbulence 
with the value of u'/U and v' jU being close to those 
of the free stream, that is 2 and 1.4 respectively, giving 
a ratio u' / v' of 1.4, close to the well known value of 

1 .2 reported for turbulent channel flow by Comtc-Bellot 
and Corrsin [8], The main point to note is the fact that 
no influence of the particles’ sizes is observed on the 
fluctuations of the velocity, but for the smallest ones 
( D < 5/m;) whose value is constant and equal to the 
value in the free stream. This fact is more clearly shown 
in Fig. 5 where profiles of the flux of particles across the 
probe volume (an indicator of the concentration of the 


Figure 7: Probality density function for o , y — 

-30; □, 0; o, 30 mm, — log-normal fit fory=-30mm. 
x=156 cm, stable stratified case. 



x (mm) 

Figure 8: Upper limit of the diffusion zone for droplets 
of size, 10 < D < 15 /tm. 

droplets) arc plotted according to their size. These pro- 
files show two very distinct regions. The first one, where 
the flux is constant, corresponds to the free stream, and 
the second one where the particle turbulent diffusion has 
already carried the droplets from the droplet-laden free 
stream to the upper part of the channel with no droplets. 
One observes that the decrease of the mass flux is almost 
exponential, and more importantly that it is, to a first or- 
der approximation, the same for all droplets whatever 
their size (note that the three straight lines are parallel 
to each other). This is further confirmed in Fig. 7 where 
we show the size PDF of the droplets measured at three 
different vertical positions. All three distributions are 
seen to be almost identical, but for a small decrease in 
the number of big droplets resulting in a small increase 
of the PDF for small diameters when y increases. 

As one can sec in Fig. 4, this was not the case at 
downstream distances closer to the splitter plate where 
a difference in PDFs is observed. Indeed, one clearly 
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Figure 9: Probality density function for o, y = 

-30; □, 0; o, 30 mm, — log-normal fit for y=30 mm. 
x=156 cm, unstable stratified case. 

notes that small droplets are diffused in the upper part 
of the channel much faster than big ones resulting in an 
important increase of the PDF for small diameters when 
y increases. The large scale eddies structure of the wake 
is responsible for this effect. The rate of diffusion of 
droplets in the upper half of the channel can be mea- 
sured from the data shown in Fig. 5. The upper limit 
can indeed be defined as the location where the data 
rale has fallen below a certain threshold. Normalizing 
these curves by the data rate in the free stream, N u , and 
using a definition similar to that used in boundary lay- 
ers, the upper limit of the diffusion zone has been cho- 
sen to be the location where the normalized data rate, 
N/N„, is 0.1. The result obtained for particles of size 
10 < D < 15 ftm, representative of the spray, is shown 
in Fig. 8. One observes that this evolution is almost lin- 
ear with a growth rate of 1%, a value very close to the 
measured normal turbulent intensity, v' /U, of 1.4%. 
Unstably-stratified case 

Similar measurements have been carried out where the 
injection of droplets was done in the upper half of the 
channel. In this second case the diffusion process is 
in the same direction as the gravity. Figure 9 shows 
three PDF taken across the test section. One clearly ob- 
serves that big droplets are diffused in the lower part at 
a much larger rate than small ones. This can be seen 
from the widening of the distribution when descending 
from y = 30 mm to y = -30 mm. Measurements are 
currently in progress to obtain a better understanding of 
the difference that exists between the stable and unsta- 
ble, the marked difference observed in the PDFs cannot 
be explained with the simple linear addition of the set- 
tling velocity of the droplets which is less than 0.1 m/s 
in all cases. Additionnal similar measurements are un- 


derway with a vertical splitter plate in order to separate 
diffusion effects from gravity effects. 
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